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We consider a two component Bose-Einstein condensate in two spatially localized modes of a 
double well potential, with periodic modulation of the tunnel coupling between the two modes. We 
treat the driven quantum field using a two mode expansion and define the quantum dynamics in 
(«~^ I terms of the Floquet Operator for the time periodic Hamiltonian of the system. It has been shown 

. that the corresponding semiclassical mean-field dynamics can exhibit regions of regular and chaotic 

motion. We show here that the quantum dynamics can exhibit dynamical tunneling between regions 
of regular motion, centered on fixed points (resonances) of the semiclassical dynamics. 
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I. INTRODUCTION 

Since the first experimental realization of a Bose condensed gas of atoms in 1995 ^ there has been considerable 
interest in the properties of Bose-Einstein Condensates (EEC's). Here we are interested in the properties of a BEG 
confined to the two spatially localized regions of a double well potential, and subjected to periodic modulation in such 
a way that the tunnel coupling between the two wells is periodically modulated. In this way it has been shown that 
the semiclassical mean- field dynamics can exhibit regions of chaotic motion . 

The primary goal of this paper is to develop an understanding of the quantum dynamics, beyond mean- field theory, 
. — of a dynamical system in a classically chaotic regime. In particular we seek quantum features of the dynamics of 
^ strongly driven BEG's that are not obtained in mean- field theory. Over the last two decades there have been numerous 
studies of nonlinear systems with a finite number of degrees of freedom; a topic that is often referred to as 'quantum 
^ chaos', see for example - In this paper however we attempt to address the quantum dynamics of a driven quantum 
(— ( ' field. In the mean- field limit, an effective classical field (a system of an infinite number of degrees of freedom) is used 
Q I to describe the condensate. The response of the system to external periodic forcing would then be described by the 
O • driven Gross-Pitaevskii equation. In the hydrodynamic limit the mean field description would look very much like a 
^ '■ forced nonlinear fluid dynamics model |^ . A similar system of equations also arise in periodically modulated nonlinear 
• i-H , optical fibers |0] . Even without external time dependent forcing nonlinear field equations of the kind considered here 
' exhibit dynamical instabilities leading to a host of interesting phenomenon for both the integrable case (e.g. solitons) 
and the non-integrable case (e.g. turbulence) |^,^. The response of condensates to forcing can encompass a large 
variety of physical situations including vortex formation in response to stirring ||lo[] , and rotating condensates in 
an anisotropic trap [Tl| ]. Adhikari Jl^ recently conducted a numerical study, in mean field theory, of two couple 
condensates with modulation of the trap frequency and also modulation of the nonlinear term due to hard sphere 
collisions. Gardiner et.al. |jl^ considered a BEG in a modulated periodic potential. 

We will consider two condensates with a tunnel coupling periodically modulated in time. In a fully quantum 
treatment the dynamics is given by a nonlinear equation for the field operator. We do not know in general how 
to deal with such systems. In this paper we adopt the usual method of a finite mode expansion, which effectively 
reduces both the classical and quantum systems to systems with few degrees of freedom. Even at this crude level 
of approximation however we can find fundamental differences between the mean-field predictions and those of the 
quantum description. 

In this paper we focus on two tunnel-coupled BEG's in a double well potential, including the nonlinear self energy 
term, with periodic modulation of the tunnel coupling. A similar model has recently been studied by AbduUaev 
and Kraenkel in the mean- field limit. They also included modulation of the energy difference between the two 
condensates as well as dissipation. They showed that, in mean- field theory, the system could exhibit chaotic oscil- 
lations of the relative population difference between the two wells. Recently, Elyutin and Rogovenko have treated a 
periodically driven double well model via the Gross-Pitaevskii equation ||2^. There are many studies of a BEG in 
a time independent double well potential p^-|20|. We will use the model of reference [Esi which is valid for small 
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condensates (few atoms, i.e. N < 1000), where one expects quantum departures from mean- field theory to be more 
significant. We make a two mode approximation for the double well system [jisl , ^ , ^ , with weak interwell tunneling. 
As particle number is a constant of the motion at zero temperature, we can use an approach based on the two mode 
bosonic realization of the SU(2) algebra. We use the Schrodinger picture and an angular momentum model for the 
quantum system in SU(2) dynamics. We find a parameter regime in which an initial atomic coherent state tunnels 
from one fixed point (resonance) to another. This behavior is analyzed from a Floquet state perspective and reveals 
an interesting feature of parity of the Floquet operator eigenstates. 

The many body Hamiltonian is specified in section ||, and our two mode approximation to it is given. Section III 
outlines the models and the corresponding equations of motion for the two mode system. The semiclassical and 
quantum results are discussed in sections [V and M respectively, followed by our concluding remarks in section VI. 



II. THE HAMILTONIAN 



The many-body Hamiltonian for an atomic BEC confined in a potential V{r, t) is given by 



H{t) 



+- Z J T ' 

— VV"^ • Vi/) -I- V)V(r, t)ij + -^V;ti/,ty3^ 
2m 2 



(1) 



where m is the boson mass, m measures the strength of the two body interaction and a is the s-wave 

scattering length, t/;^ — 'i/j'^{r,t) and ip = 'ip{r,t) are the Heisenberg creation and annihilation operators of particles 
at position r. V{r,t) — Vx{r){l + ecoswut) -I- Vy{r) + Vz{r), and e < 1. V{r,t) is a time dependent potential, 
parameterised by the modulation strength e, and the driving frequency lojj. Throughout this paper e = implies 
non-driven dynamics and e > corresponds to driven (modulated) dynamics. 

We will follow the model of reference in which the trap potential is taken to be a symmetric double well in the 
x-direction and harmonic in the y, z-directions. 
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(2) 



where LOt is the trap frequency in the y—z plane, and where h gives the strength of the confinement in the x-direction. 
The potential has stable fixed points at ri = +'i'oX, V2 = —qo^ near which the linearised motion is harmonic 
with frequency t^o = 'Zo(86/m)^/^. For simplicity we will set uJt — ijJq- The length scale is set by the rms position 
fluctuations, rp, in the ground state of the harmonic potential near the fixed points; = •v/?i/2mwo- In terms of this 
parameter the barrier height separating the two wells may be written B = (?icJo/8)(<?o/''o) • 

We assume the potential is such that there are two nearly degenerate single particle energy eigenstates below the 
barrier and expand the field in terms of two localised single particle states, uj^{r), at each stable minima; 



Ud{r) =ito(r-rrf), d= 1,2 



(3) 



where uo(r) is defined to be the normalized single-particle ground state mode of the harmonic potential near each 
stable fixed point, r^;, each with energy Eq. The details can be found in [Q. The energy eigenstates of the double well 
may then be approximated as symmetric (-I-) and antisymmetric (— ) combinations of the localised states with energy 
eigenvalues given in first order perturbation theory by E± = Eq ± hn/2 where the energy level splitting determines 
the tunnelling frequency and is given by 



(4) 



As discussed in |18 the validity of the two mode approximation requires fi/tJo << 1- 

The effect of the modulation is to directly modulate the parameter h so that it becomes 6(1 -I- ecosojci). This has 
two important effects; it modulates the harmonic frequency around each fixed point and it modulates the tunnelling 
frequency, so that these quantities become time dependent wo(Oj ^(0; where 



n{t) — rio(l + e cosLUnt) exp 



-^Wl + ecosLJDt-1) 



(5) 
(6) 



The two mode approximation then results in the following Hamiltonian p8[ 
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H2[t) = hLUo{t){c\ci + 4c2) + ^^(Cl4 + C\C2) 

+ hK{{clfcl + iclfcl), (7) 
Cd{t) = / d^ru*a{r)Mr,t), (8) 



where 



and similarly for cj;(t), d = 1,2. cj,ci and cJ,C2 are creation and annihilation operators for atoms in the local 
modes in wells 1 and 2 respectively. The remaining parameters in the Hamiltonian are as follows, k represents the 
particle-particle interaction and measures the tunnelling rate of atoms between wells 1 and 2. We will assume 
that the system remains in a total number eigenstate, with eigenvalue N. We first note that we will be interested in 
modulation frequencies that are of the order of the unmodulated tunnelling frequency, lud ~ ^o- The condition for 
the validity of the two mode approximation then ensures that lud << ujq. We can then replace the first term in the 
Hamiltonian by a slowly varying number hu!o(t)N, which has no influence on the quantum dynamics. The temporal 
modulation of the second term, which is responsible for tunnelling, has a strong influence on the dynamics as we show. 
If e is small enough we can approximate the modulated tunnelling frequency by 

i}{t) fa ilo{l + ecosLUot) exp[— je cos uj at] (9) 

where 7 = go/(4^o)- While e may be small we cannot necessarily claim that 76 is small. We may expand the 
approximate expression for Q{t) using Fourier series, 

CXD 

n{t) fa iloil + ecosLUDt) J2 In{-Ji)e-""^''* (10) 



where /„ (x) is a Bessel function. In general we see that the direct modulation of the potential leads to a non-harmonic 
modulation of the tunnelling frequency. In what follows we will neglect the higher harmonics of the modulation 
frequency, and consider the time dependence of Q{t) to be simply 

n{t) ^^n^l + e cos UDt) (11) 

While we have assumed that e is small, this constraint is not as strong as one might expect. When we study the 
ratio of the approximate modulated tunnelling to the complete expression given in Eq.(|6|), 0,{t)/il{t), we find that 
for appropriate values of Qo/rQ, the ratio is close to unity even for values of e as large as 0.3. Henceforward we will 
assume the simpler harmonic modulation of Q{t) — f2o(l + ecosw^i). 

The separation of the higher excited states from the nearly degenerate ground states is of the order of ojq. As the 
validity of the two mode approximation requires that Qq << loq, and we assume that the driving frequency is of 
the order of the tunnelling frequency. The driving frequency is very far from resonance with higher excited states of 
the double well and we can continue to make a two level approximation when the modulation is present. Of course 
when the tunnelling and the modulation is slow, the interesting dynamics will occur on a long time scale and one may 
legitimately ask if dissipation should not be included. In this paper however we will neglect the effect of dissipation 
and focus only on the nonlinear quantum dynamics. Our results are thus strictly only valid for zero temperature. The 
validity of an expansion in terms of localised single particle states places a restriction on the atomic number |l^ . We 
are thus only concerned with very small condensates at temperatures low enough that the non-condensate fraction, 
and the resulting dissipation, can be neglected. It is only in this regime that quantum (rather than semiclassical) 
nonlinear effects in driven EEC's could be observed. 

Using the generators for SU(2) dynamics, 

(12) 
(13) 
(14) 

the two mode Hamiltonian may be written. 
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H2{t) = hn{t)j, + 2nKj^. (15) 

It is worth noting the symmetry of this Hamiltonian under parity reflection of — > —Jx- The dynamics of this two 
mode system are greatly simphfied by this symmetry. 

The operator represents the difference in atom number between wells 2 and 1 while the total number of atoms 
is given by c\ci + c\c2 — Ni + N2 — N . As a consequence of number conservation, the operators of Eqs.(|l2|,|l3|,^4|) 
have a Casimir invariant given by 

Jy and Jz are of no particular significance for our present purposes and are not considered further. 



III. THE EQUATIONS OF MOTION 

In reference |l^ the mean-field equivalent of the two mode approximation used here was shown to result from a 
moment factorization ansatz in the Heiscnberg equations of motion. In effect this is equivalent to a classical two mode 
approximation for the Gross-Pitaevskii equation. We start by using the Heisenberg equation of motion to find the 
time evolution of the three operators in Eqs.( [l2| , |l3| , |l4[) . Using the familiar commutation relations [Ji, Jj] = —iJk, we 
obtain three equations of motion for the operators Jx, Jy and Jz- 

^ = -mJy, (17) 
^ = n{t)Jx-2K(jJx + JxJz), (18) 

^ = 2K{JyJx + JxJy). (19) 

We transform these Heisenberg equations of motion into a semiclassical system of equations by taking the expectation 
value of the above equations and factorizing operator products such as (JxJz) = {Jx){Jz)- This leaves us with classical 

numbers in place of the J operators. We make the identification Si = we scale the J equations by the total 

number of atoms N. We now write down the semiclassical system of scaled equations, 

Sx = -n{t)Sy, (20) 

Sy^nit)Sx-inNSxSz, (21) 
Sz = 4.KNSxSy, (22) 



with the usual constraint to the sphere. 



Si + Si = 1/4 = 5*2. (23) 



We identify three real variables Sk, corresponding to the operators J^, 

1 
2 



Sx = ]:{blh2-blbi) (24) 



Sy^^(blb2~b,hl) (25) 

Sz^\{blb2+bibl) (26) 

where the bdS are c-numbers, and are the expectation values of the second quantized annihilation and creation 
operators Cd- These real variables Sk can be shown to represent the following physical parameters, NSx represents the 
difference in atom number between wells 2 and 1, NSy represents the momentum of the condensate and N Sz represents 
the population difference between the global symmetric and antisymmetric modes of the confining potential pa . 
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IV. SEMICLASSICAL DYNAMICS 



While wc can imagine a vector tracing out some trajectory on the surface of the sphere, see Fig. |l|., it is convenient 
to make a stereographic projection of the dynamics onto a two dimensional plane. 

We use a stereographic projection of a sphere onto the plane z = 1/2 via the transformation, 



iSqi 



1/2 + 

Using the definitions a — Reiy) and b = Im{v) we write, 



a = 



1/2 + 



(27) 



(28) 



1/2 + 5"^ 

Setting A = nN/flo, we rewrite the system of Eqs.( pO| , pl 22) as, 

a = -(1 + €cos-^t)b - 4A -— 

no l + aM 

6 = (1 + ecos fracuj DnQt)a — 2Xa 



(29) 



fe2 



1 



62 



(30) 



Where we have scaled the time t by ilo and we set ud = 1.37f2o- The equator of the sphere is mapped to the unit 
circle in (a,b) phase space due to our choice of radius (r = ^) for the sphere and the position of the wells are mapped 
to (±1,0). 

For e = 0, i.e. ri(t) = Hq, the system (30) is integrable. If A < 1/2, the only feature is a center at the origin. At 
A = 1/2, the origin undergoes a pitchfork bifurcation creating two stable centers, which for A > 1/2 have position 



2A- 1 

2A + r 



5 = 0, 



and are surrounded by a 'figure 8' type separatrix. The integral curves for this system are. 



c(a, b) 



(1 



2\2 



-52) 



(l + a2 + 52) + 2Aa2 



(31) 



(32) 



On the separatrix c = 1, for c < 1 solutions to Eq.(|3^ lie inside the separatrix and for c > 1 they lie outside. 

Fig. Ij. shows the phase plane for e = 0, illustrating the integral curves. If we consider starting in one well, say 
(a = +1,5 = 0), then for A < 1, complete population oscillations between the wells occur |l^. At A = 1 we find 
the separatrix intersects the centers of the wells and complete oscillations no longer occur. As A increases beyond 1 
(larger N, i.e. more atoms) the solutions of Eqs.(|3^) become localized near the fixed points of the phase plane. This 
is known as self-trapping. Physically, this can be interpreted as an increase in the self energy with particle number 
which modifies the effective single particle potential of the system and reduces the single particle tunneling |l^ . 

A very similar model with harmonic driving is considered by Elyutin and Rogovenko [ pT[ . In that paper they 
find fixed points and a separatrix for their model analagous to the fixed points and separatrix in the model treated 
here. They consider the break down of the regular phase space due to the overlap of resonances, and show tunnelling 
between the classical fixed points of the semiclassical phase space. We go a step further here in that direction firstly 
to analyze the dynamics of fixed points of the Poincare map (resonances) and initial conditions localized on those 
semiclassical resonances and secondly to analyze a full quantum description. 



Resonances 



For e > the modulated tunneling rate n(t) gives rise to resonances. Outside the separatrix, a resonance will 
appear when there is a rational ratio between the driving period and the period of motion of the system on a given 
integral curve. We know on a particular integral curve c(a, 5), the system period is given by. 
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AK{k) 



where K{k) is a complete elliptic integral of the first kind, — 



(33) 



W^—W- 



and , w_ are given by 



w± = 



-(2A - 1) ± V(2A- 1)2 +8A(c- l)/c 



(34) 



In particular, the integral curve given by Tsys{c,X) = is resonant, breaking up under perturbation to give two 
stable period- 1 resonances which for i = lie on the momentum axis. These can be seen in the Poincare section. 



t = mod{ " ), which here is simply a stroboscopic map taken at i„ 



, 2irao 



In such a stroboscopic map, see 



Fig. 1^., the resonances appear as fixed points. They lie outside the separatrix, which has been replaced by a 'chaotic 
sea', and are distinct from the two period-1 fixed points nearer the origin. For ^ > 1, here it was taken as 1.37, 
these period-1 resonances exist for A > (^)^ — l)/2. 

In order to make a comparison with the quantum dynamics we need to consider an initial classical distribution 
rather than a single point, and compute average values of the dynamical variables as the distribution evolves in time. 
The classical phase space density (distribution) we require to simulate the quantum state must correspond to a system 
of points localized on the spherical classical phase space and should describe the optimal simultaneous determination 
of the phase space variables, that is the components of angular momentum. These states are known as the SU(2) 
coherent states ||2^ and are defined in Eq.(|^). 

For a quantum state \4>), the distribution that describes the output for optimal simultaneous measurement of 
angular momentum is 



where a = e 
is M 



P(a) = |(«|^)P (35) 
"^tan(|). If the system begins in an initial SU(2) coherent state |ao): the corresponding distribution 



P(a|ao) 



(1 -I- aoa*)(l -I- a^a) 
(l + |aoP)(l + |aP) 



2i 



(36) 



where j is the total angular momentum quantum number. This is the classical distribution that we sample to compute 
the ensemble dynamics. The distribution is normalized with respect to the spherical integration measure, which in 
the stereographic coordinates is d^Q-^i _(_ |q,|2~)-i_ if ^e start an initial distribution of points in a stable resonance 
island, the population remains confined by the invariant tori bounding the island, see Fig. |^. In the next section we 
shall see how the equivalent quantum simulations produce a very different result. 



V. QUANTUM DYNAMICS 



It is difficult to solve nonlinear Heisenberg equations of motion so we work in the Schrodinger picture using an 
orthonormal basis defined by the eigenstates of J^, as in [^. Using the familiar angular momentum notation we 
write our state \ip) as, 

1^^)= ^ c„(<)b-,n) (37) 

n=-j 

and substitute into the Schrodinger equation to get 

j 

ihcn{t) = {j,n\H\j,m)cm{t)- (38) 

This is a set of linear coupled equations (2j -I- 1 of them) and we use a Runge-Kutta numerical routine to solve them 
for the coefficients c„(i). Given these coefficients we can obtain the various expectation values (Ji) and thus the 
solution of our model. We once again treat the driven and undriven systems separately. For e = we obtain regular 
collapse and revival sequences as shown previously in fl^]. 
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For e > the system is described by a time periodic Hamiltonian. In that case the appropriate description is in 
terms of the Floquet operator F, of the system which maps the state from one time to a time exactly one modulation 
period later 

^ „ = 0, 1, 2, . . . . (39) 



The Floquet map is the quantum equivalent of the classical Poincare section defined in section [V. We obtain the 
Floquet Operator in the basis which diagonalizes Jz as follows. Solve the Schrodinger equation over one modulation 
period of time T to obtain \ijjrniT)) for the set of initial conditions to) to = —j, —j + 1, . . . , j — 1, j} . The solution 
vectors, ItpmiT)), form the columns of the Floquet Operator matrix F. 

In Fig. ||., we plot the mean value of for an initial state localized on the same period-one resonance as in Fig. IJ. 
We use the SU(2) coherent states as the initial localized states. These are defined on the sphere by [fsl, 



E 



^ ^ (40) 



m + j J (l + |a|2)J 



where a — e*'^tan(|), and 9 and (f> are the spherical polar coordinates. For the parameters e = 0.3 A — 1.6 it is 
found that after 310 modulation periods, the system state appears to be localized on the other resonance, see Fig. ||. 
This behavior is known as dynamical tunneling . Dynamical tunneling is a different process from the well-known 
case of barrier tunneling. There is no energetic or potential barrier. None the less in dynamical tunneling the system 
exhibits classically forbidden motion, motion which cannot be accessed by a classical particle moving on any trajectory. 
Dynamical tunneling in a single atom system has recently been observed by two groups p8| , ^ . 

We can understand the tunneling in the following way. The Hamiltonian is invariant under rotations about the 
z-axis, i.e. operations of e'^"'^ leave the Hamiltonian unchanged. Thus the Floquet states, by the nature of their 
construction must fall into two sets, corresponding to the eigenvalues of e™'^' . These two sets contain vectors of odd 
and even parity. 

Let be a state localized on one of the period-one resonances with coordinates (—1.72, 0), while \tp+), is localized 
on the other period one resonance, (-fl.72, 0). These states are related by a rotation of \tlj-) by ±7r about the z axis, 

e"-^^|V'+) = |V-). (41) 

As the Hamiltonian is invariant under rotations of tt around the z axis, the eigenstates of the Floquet operator 
fall into two parity classes defined by this rotation. If the initial state is well localized on one of the period-one 
resonances, it will tend to have maximal support in a two dimensional subspace spanned by two particular Floquet 
eigenstates of opposite parity. We label these two simultaneous eigenstates of F and parity as \4>±), where e"^± are 
the corresponding eigenvalues of the Floquet operator. We then write the 'negative' and 'positive' localized states, 
\ip±) as a superposition of these parity eigenstates states, i.e. 

= fci|0+)+fc2|<^-), (42) 

and, 

I =fci !(/)+) -fc2|0-), (43) 

where /ci, ^2 are very nearly equal and together almost exhaust the normalization condition, /c^ -f fc^ « 1. 

For complete tunneling we must have, \tp+) — > How many iterations of Eq. (|39|) satisfy this tunneling condition? 

Suppose we start in the state After n iterations of the Floquet map we reach the state 



1^") = F"|V^+) 



-^(e-^+|0+)-He-^i0_)) 



(44) 



When the relative phase of the two states in this superposition is shifted by tt, the state becomes almost equal to a 
state localized in the other period-one resonance lip-)- This determines the tunneling time as 

0- - 0+ 

To find the particular Floquet eigenstates on which an initially localized state has maximum support we simply expand 
the initial state in the basis of Floquet eigenstates. This also gives the corresponding eigenphases and thus the 
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tunneling time. By this method we find the tunneling time should be nt = 310. This compares favorably with the 
tunneling time of about n = 313 evident in Fig. |^. 

We expect to see the variance of Jx at a minimum whenever the component is most localized around each of 
the resonances, which is when the state should be close to a minimum uncertainty state. At the half way point, a 
coherent superposition of two states of equal and opposite mean values of should occur, and the variance should 
be a maximum. This behavior is evident in Fig. ^ 



VI. CONCLUSION 



In this paper we have used a two mode approximation to describe the quantum dynamics of a BEC in a double 
well potential with periodic modulation of the potential barrier separating the two wells. In mean-field theory (Gross- 
Pitaevskii limit) we find, for appropriate parameter values, a mixed phase space of regular stable motion, associated 
with fixed points, coexisting with chaotic dynamics. We have shown the existence of dynamical quantum tunneling 
between period-one fixed points (resonances). This dynamical tunneling is distinct from the single particle tunneling 
present in the undriven system and may be interpreted as a new signature of the quantum dynamics of a nonlinear 
quantum field. 

Double well potentials similar to those modeled here have been experimentally realized, see for example |l|,||,|§. 
While current experiments are unlikely to be well described by the two mode model given here, we expect dynamical 
tunneling to be sufficiently robust that an experimental realization of the dynamical quantum tunneling predicted here 
can be observed. The recent observation by the NIST group of dynamical tunneling of a single particle system 
in a modulated periodic potential in fact used a BEC, although not in a regime where the mean-field is significant. It 
would not be difficult to modify the experiment so that the nonlinear effect of the BEC would be significant. 
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FIG. 1. Diagram illustrating construction of the stereographic projection from a sphere. Points 1 and 2 represent the 
positions of the two wells that constitute the double well potential we are modeling, p is a point on the sphere projected to the 
point P on the plane, a and b are the coordinates in the plane, as defined in the text. 
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FIG. 2. Phase space plot of the semiclassical dynamics showing period-1 fixed points and separatrix. e = 0.0, A = 1.6. 
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FIG. 3. Phase space plot of semiclassical driven dynamics showing period-1 resonances at o = ±1.72, 0.0. e = 0.3, A = 1.6, 
WD = 1.37. 
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FIG. 4. Plot of the weighted mean of semiclassical population difference < Sx > and variance Var < Sx >, with initial 
distribution centered on the resonance a = 1.72, b = 0.0. e = 0.3, A = 1.6, ujd = 1.37. 
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